# -*- coding: utf-8 -*-
###################################################
#    This file is part of blockIMH.
#
#    blockIMH is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    blockIMH is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with blockIMH.  If not, see <http://www.gnu.org/licenses/>.
###################################################
#! /usr/bin/env python

from __future__ import division
import os
import sys
import scipy.weave as weave
from scipy.stats import norm, cauchy, uniform
from scipy import special
from numpy import random, exp, sqrt, power, log, newaxis, zeros, \
        zeros_like, float32, int32, repeat, minimum, ones, mean, \
        average, sum, var
## In fact CUDA is not used here
try:
    import pycuda
    import pycuda.autoinit
    import pycuda.gpuarray as gpuarray
    import pycuda.driver as cuda
    import pycuda.cumath as cumath
    from pycuda.compiler import SourceModule
    CUDA_ENABLED = True
except:
    CUDA_ENABLED = False
from IMH import blockIMH, simpleIMH
import permutations
from localfolder import get_path
#import rpy2.robjects as robjects

## Function of interest
## (here we take the identity, so 
## that our estimators will estimate
## the mean)
def h(x):
    return (x)

#####################################################################################################
## 

## Nb of parallel processors
try:
    for arg in sys.argv:
        arg = arg.replace("\s+", "")
        if "p=" in arg:
            p = int(arg.replace("--p=", ""))
        if "b=" in arg:
            b = int(arg.replace("--b=", ""))
        if "nbsim=" in arg:
            nbsim = int(arg.replace("--nbsim=", ""))
except:
    raise NameError("""
problem processing the arguments: the only available options are:
--p=X where X is an integer denoting a block's length,
--b=X where X is an integer denoting the number of blocks (hence T = b * p is the chain's length),
--nbsim=X where X is an integer denoting the number of independent runs to perform.
(don't put spaces around the "=" signs)
""")
if not vars().has_key("p"):
    print "no p given in argument, default is p = 4"
    p = 4
if not vars().has_key("b"):
    print "no b given in argument, default is b = 1"
    b = 1
if not vars().has_key("nbsim"):
    print "no nbsim given in argument, default is nbsim = 100"
    nbsim = 100
## Nb of iterations
T = p * b
x_0 = 0
## Nb of independent runs

print """Starting block IMH algorithm with:
    p = %i
    b = %i
    T = %i
    %i independent runs
""" % (p, b, T, nbsim)

#####################################################################################################
## 
## Where are we?
basefolder = get_path()
## Set some filenames to save the results
resultsfolder = os.path.join(basefolder, "results")
# Allocate memory to save the results
results = {}
results["basicestimvector"] = zeros(nbsim)
results["blockestimvectorCircPermut"] = zeros(nbsim)
results["estimRB1vectorCircPermut"] = zeros(nbsim)
results["blockestimvectorNoPermut"] = zeros(nbsim)
results["estimRB1vectorNoPermut"] = zeros(nbsim)
results["blockestimvectorRandPermut"] = zeros(nbsim)
results["estimRB1vectorRandPermut"] = zeros(nbsim)
results["estimRB2vectorRandPermut"] = zeros(nbsim)
results["blockestimvectorHalfRandPermut"] = zeros(nbsim)
results["blockestimvectorStratRandPermut"] = zeros(nbsim)
results["blockestimvectorTempPermut"] = zeros(nbsim)
results["estimIS"] = zeros(nbsim)

#####################################################################################################
## Main part

# Profile the computer intensive part
import cProfile
cProfile.run("""
for simu in range(nbsim):
    print simu
    imh = blockIMH(x_0, T, p, h)
    results["basicestimvector"][simu] = imh["basicestim"]
    results["blockestimvectorCircPermut"][simu] = imh["blockestimCircPermut"]
    results["estimRB1vectorCircPermut"][simu] = imh["estimRB1CircPermut"]
    results["blockestimvectorNoPermut"][simu] = imh["blockestimNoPermut"]
    results["estimRB1vectorNoPermut"][simu] = imh["estimRB1NoPermut"]
    results["blockestimvectorRandPermut"][simu] = imh["blockestimRandPermut"]
    results["estimRB1vectorRandPermut"][simu] = imh["estimRB1RandPermut"]
    results["estimRB2vectorRandPermut"][simu] = imh["estimRB2RandPermut"]
    results["blockestimvectorHalfRandPermut"][simu] = imh["blockestimHalfRandPermut"]
    results["blockestimvectorStratRandPermut"][simu] = imh["blockestimStratRandPermut"]
    results["blockestimvectorTempPermut"][simu] = imh["blockestimTempPermut"]
    results["estimIS"][simu] = imh["estimIS"]
""", "prof")

print "Finished!"
print "Now some profiling..."
# Print profile statistics
import pstats
pS = pstats.Stats('prof')
pS.sort_stats("cumulative").print_stats(10)
pS.sort_stats("time").print_stats(10)
pS.sort_stats("cumulative").print_stats("omega")

import cPickle
picklefilename = "results/results-p%iT%in%i.cpickle" % (p, T, nbsim)
f = open(picklefilename, "w")
cPickle.dump(results, f)
f.close()
import cPickle
import pickletoRdata
pickletoRdata.pickle2RData(picklefilename)
os.remove(picklefilename)

print "Goodbye!"

